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UNSTEADY  TRANSONIC  FLOWS:  TIME- LINEARIZED  CALCULATIONS 


A.  Richard  Saabaas  and  K.-Y.  Fung 

Aerospace  and  Mechanical  Engineering,  University  of  Arizona,  Tucson,  Arizona 


An  accurate  and  efficient  method  of  computing  unsteady  transonic  flow  is 
described.  The  flov  is  Linearized  about  an  experimentally  measured  or  numeric¬ 
ally  calculated  steady  state,  as  represented  by  a  given  pressure  distribution. 

For  a  given  mode  of  motion,  the  amplitudes  and  phase  lags  of  the  lift  and 
moment  coefficients  at  a  given  reduced  frequency  are  found  by  superposition 
from  an  indicial  response.  The  computational  effort  is  reduced  by  treating  shock 
waves  as  discontinuities,  and  by  applying  the  correct  linear  far-field  behavior. 

A  novel  method  of  modeling  the  Indicial  response  provides  an  analytical  formulas 
for  the  dependence  of  Che  amplitude  and  phase  lag  on  the  reduced  frequency. 


INTRODUCTION 

A  combination  of  technical  advances  should 
improve  the  fuel  efficiency  of  transport  aircraft 
by  fifty  percent  In  the  next  decade.  Analogous 
improvements  in  the  transonic  performance  of  mili¬ 
tary  aircraft  should  also  be  realized.  These  large 
gains  will  come  from  a  combination  of  improvements 
in  engine,  structural,  and  aerodynamic  efficiency. 
More  than  half  will  come  from  Improvements  In  the 
aerodynamic  efficiency,  including  active  control, 
and  the  use  of  composite  materials  in  the  primary 
structure.  Part  of  the  Improvement  In  aerodynamic 
efficiency  will  result  from  flight  at  supercritical 
Mach  numbers  with  subcrltlcal  levels  of  lift  to 
drag  ratio  and  high  lift  coefficients  at  near  sonic 
flight  conditions. 

This  improved  transonic  performance  mandates 
an  accurate  prediction  of  aeroelastlc  behavior  at 
transonic  Mach  numbers.  Cf  special  concern  are 
flutter  boundaries.  In  1976  Farmer  and  Hanson  (1) 
reported  chat  the  flutter  boundaries  of  two  dynam¬ 
ically  identical  wings  were  markedly  different  at 
transonic  Mach  numbers  due  to  very  minor  differ¬ 
ences  in  wing  profile  chickness.  The  results  of 
their  measurements  are  shown  in  Figure  l,  indica¬ 
ting  the  reduced  flutter  boundary  for  the  wing  with 
a  "supercritical"  profile. 

Today  we  understand  well  the  qualitative 
behavior  of  invlscld  steady  and  unsteady  transonic 
flows,  and  we  have  rudimentary  understanding  of 
viscous  effaces.  For  flight  regimes  that  Involve 
unseparaeed  flows  tha  main  lngredlanc  in  the  calcu¬ 
lation  of  fluttar  boundaries  Is  an  accurate  deter¬ 
mination  of  the  steady  pressure  field.  This  may  be 
determined  either  by  experiment  or  by  calculation. 
But  once  it  la  known,  the  response  of  the  vlng  to 
pitching,  plunging  or  aileron  motion  may  be  found 
by  numerical  means.  Tha  most  essential  lngredlanc 
In  predicting  this  response  Is  an  accurate  predic¬ 
tion  of  the  motion  of  any  shock  waves  present  in 
the  flov  (see,  e.g. ,  (2)).  Numerical  algorithms 
chat  capture  shock  waves  must  use  relatively  fine 
grid  spacing  near  the  shock  wave  if  they  are  to 
predict  Its  motion  due  to  the  smell  changes  of 
lneereet  In  flutter  studies.  But  this  shock 
motion  may  be  predicted  accurately  by  a  time- linear¬ 
ized  algorithm  using  a  relatively  coarse  grid  If 
the  calculations  are  done  correctly.  This  has  not 
generally  been  the  case,  with  ocher  Investigators 
ignoring  this  essential  efface  (3,4,3). 


tfa  report  here  on  our  two-dimensional,  cime- 
linearlzed  computations,  which  not  only  properly 
account  for  shock  wave  motion,  but  are  able  to 
resolve  them  even  though  the  grid  used  to  calculate 
the  flow  is  relatively  coarse.  In  order  to  avoid 
the  reflection  of  the  unsteady  disturbances  from  the 
grid  system,  only  a  moderate  amount  of  grid  stretch¬ 
ing  Is  employed  away  from  the  airfoil.  To  avoid 
unnecessarily  large  computational  domains,  the 
linearized  far-field  for  an  unsteady  vortex  with  a 
circulation  determined  by  the  airfoil's  lift  is  used 
to  evaluate  the  potential  there  (6).  The  airfoil's 
response  to  a  given  mode  of  motion  is  determined  by 
superposition  from  chat  for  an  indicial  motion.  In 
many  cases  this  indicial  response  can  be  modeled  In 
a  simple  way,  providing  an  analytic  result  for  the 
dependence  of  the  lift  and  moment  coefficient's 
amplitudes  and  phase  lags  on  reduced  frequency.  A 
novel  feature  of  this  modeling  is  that  of  a  sequence 
of  harmonic  oscillators,  each  of  which  improves  the 
previous  simulation  of  the  indicial  response.  This 
provides  an  analytical  formula  for  tha  dependence 
of  the  amplitudes  and  phase  lags  on  the  reduced 
frequency  (7). 

The  computational  efficiency  of  the  time- 
linearized  calculation  of  the  amplitudes  and  phase 
lags  for  a  range  of  reduced  frequencies  Is  compared 
with  nonlinear  and  frequency  domain  computations. 

Tha  time- linearized  calculation  of  an  indicial 
response  can  result  la  a  factor  of  ten  or  more 
reduction  in  computational  effort;  this  is  especially 
significant  in  three  dimensions.  The  modeling  of 
the  indicial  response  by  a  sequence  of  two  harmonic 
oscillators  can  also  reduce  computational  effort. 

GOVERNING  EQUATIONS 

As  noted  above,  time- linearization  about  a 
known  steady  state  Is  an  effective  mechanism  for 
datsrmlnlng  the  unsteady  responsa  to  a  given  mode  of 
motion.  For  flows  chat  are  unseparaeed,  and  not  at 
Incipient  separation,  an  invlscld  treatment  of  the 
unsteady  flow  should  be  adequate  for  flutter  studies. 
It  la  li^ortaat,  however,  that  the  nonlinear  and 
viscous  aspects  of  the  underlying  steady  state  be 
determined  accurately.  This  may  either  be  done  by 
experiment,  or  by  a  reliable  computational  algorithm 
such  as  Cr unfoil  (8).  In  any  case,  ve  assume  chat 
an  accurate  steady  state  pressure  distribution  has 
been  used  to  provide  the  input  fer  the  inverse  cal¬ 
culation  of  the  airfoil  shape  that  will  provide  this 
pressure  distribution  when  the  steady  stats  flow  is 


computed  using  eh«  small  perturbation  approximation. 

As  tin  ce  al.  (9)  obsarvad  more  than  thirty 
years  ago,  within  tha  contaxt  of  tha  small  parturba- 
tion  approximation  tha  basic  aquation  govaming  the 
unsteady  motions  is  linaar  unless  the  reduced 

frequency,  k,  is  0(6  ).  Here  5  is  tha  maasura  of 

o  o 

the  perturbation  potential  and  k  •  wc/tT,  where  cu  is 
the  frequency,  c  the  chord  and  U  the  freestreaa 
speed.  Thus,  if  we  write  the  velocity  potential, 

♦  .  as 

<•  ■  Uc{x  4*  5 

o 

where  }  is  the  perturbation  potential  and  5  Is 

o 

some  measure  of  the  disturbance,  the  governing 
aquation  is 

-  (2k/«9)t^*jte  +  {(1  -  tr>/«0  -  (r  +  DM2^}*^  (1) 


2 

where  terms  G(k  4  . /5  )  and  Q(kd  p  )  have  been 

tt  O  C  XX 

neglected  because  k  *  0(5  ).  (Hare  the  time  has 
o 

been  nondimenslonalizad  by  tha  circular  frequency 
and  the  spatial  coordinates  by  the  airfoil  chord, 
and  the  y  coordinate  is  scaled  by  $o1/2.)  The 

second  is  of  little  consequence.  Neglecting  the 
first  is  equivalent  to  disregarding  one  of  the 
characteristics  and  assuming  disturbances  propagate 
downstream  at  infinite  speed,  and  has  important 
computational  advantages.  The  boundary  condition 
at  the  body  Is 

3/2 

50  *  «'{*)  *  3V*-e>  ♦  MT  (*,e)  (2) 

where  the  body  is  given  by 

y  -  5Y(x)  *  5T(x,c>. 

Across  the  airfoil  wake  the  Jump  in  the  pressure 
coefficient  must  vanish.  Thus, 

([  H  C?(x,0,c)]l»  -  £px(x*0,c) 

+  k^c(x,o,t)]I  -  0, 

where  j[  (■••)][  indicates  the  jump  across  the  wake. 
This  implies  that  in  the  wake 

dd(x,o,c)  -  r(x  -  kt)/Uc5  .  (3) 

o 

In  both  these  boundary  conditions  we  have  retained 
terms  of  0(k) ,  which  is  not  consistent  with  the 
approximation  made  in  Eq.  (1).  But  as  Refs.  (10) 
and  (11)  demonstrate,  this  gives  good  agreement 
with  the  results  of  linear  theory  for  values  of  k 
up  to,  and  evan  abova,  1.0.  This  ttrm  is,  of 
course,  also  retained  in  tha  evaluation  of  the 
pressure  coefficients.  As  noted  in  Ref.  (2),  tha 
appropriate  measure  of  5  is  max(53^2, 5^/2  (k«>3/2) 
3/2  0 

and  normally  5  .  Far  from  tha  airfoil  (see  (6)), 

♦Ot.y.t-)  -  -  e;>^gi2.  d«;  (4) 

where 

f(x,yf,tf)  *  H(t*  +  x  -  /(x2  *  y,2))x 

.'mb’1  *  2xf'  Z^lh  t  si 

r 


-  -1  v' fe'2  +  2xt  -v’2)  -  t’ , 

•  can  y*  - >  * 

and 

c’  -  tU  -  M2) /icM2  ,  r  -  y/Cl  -  M2)/^. 

Here  H  is  the  Heaviside  unit  step  function. 

£n  addition  to  Eq.  (1)  and  the  boundary  condi- 
tlons  (2)-(4),  a  shock  jump  condition  needs  to  be 
Imposed  if  the  shock  wave  Is  to  be  treated  as  a 
discontinuity  rather  than  "captured"  by  the  numeri¬ 
cal  calculations.  Because  the  former  is  tha  Intent 
here,  we  need  to  note  that 

-  (2kM2/4a)j[*x]f(d*/dt)j 

-  Cl  -  *£>/«,  -  <Y  ♦  UM2*^  IE*  jf  +I*732  •  0,(5.) 

on 

(dy/dx)#  *  -  (5b) 

where  [£  (•••)][  and  (*  *  * )  indicate  the  jump  in  and 
average  of  (•••)  across  the  shock  wave.  Equation (5a) 
insures  the  conservation  of  masa.  The  conservation 
of  momentum  is  replaced  by  the  lrrotationallcy  con¬ 
dition  Eq.  (5b)  or  its  equivalent,  [£$][■  0. 

Time- linearization 

We  use  the  ADI  technique  Introduced  by  BaJIhaus 
and  Steger  (12)  to  compute  the  steady  state  solution 
of  Eq.  (1),  so(x,y),  subject  to  the  steady  boundary 

conditions  Implied  by  Eqs.  (l)-(4),  using  the 
coordinate  stretching  of  Ref.  (13).  Aside  from  the 
far-fleld  condition  (4)  and  the  inclusion  of  terms 
of  0(k)  relative  to  0(1)  in  (2)  and  (3),  this  is 
equivalent  to  NASA  Ames  computer  code  LTRAN2  of 
Ballhaus  and  Coorjian  (14).  We  next  linearize  about 
this  steady  state  by  assuming  that 

»(x,y,t)  «  4  (x,y)  +  (4/43/2)i(x,y,c)  *o(«/«3/2), 

OO  0 

*  3/2 

where  5/5q  *  o(l).  This  gives,  with  5  ■  5  , 

-  (2k«f/«0)ixe  ♦  (Qi  - 

-  « *  %>* *  v  •  °*  *•> 

with 

4y(x,o,t)  -  Y^x.t)*  kYt(x,t)  (6b) 

and 

!E*x<X,o,t)  +  kdc(x,o,t) JJ  -  0.  (6c) 

With  the  linearization  of  the  solution  about  a 
steady  state  at  time  t  *  0,  we  use  (4)  for  the 
potential  far  from  the  airfoil  as  the  circulation 
departs  from  its  steady  state  value. 

As  noted  earlier,  tha  proper  account  shock 
motions  ars  of  prime  importance  in  unsteady  transo¬ 
nic  flow.  We  thus  use  the  procedure  of  Ref.  (13) 
to  account  for  shock  motions.  Because  the  shock 
waves  are  nearly  normal  to  the  frees cream  we  assume 
that  this  Is  tha  case  and  satisfy 

!♦  3  -  o 


2 


on  the  normal  shock  approximation  to  Eq.  (5),  viz., 

L*_k  ,  + ; }. 

2k  (Y  +  DM2  ** 


dt 


Again  va  Un*ari2e  about  the  steady  state,  writing 
In  this  approximation 

*  <o  *  * +  ($/$y2)x(t> 

S  09  0 


which  gives 

dX(C)  .  Y£  1  ♦x(*<w.o.c)  (7) 

as  the  equation  that  keeps  track  of  the  shock  wave 

Sos it Ion.  Straightforward  linearization  of 

r  d  +■  (5/6*' 2)  ill  *  0  gives  the  expression  chat 
“0  o  •** 

determines  }  behind  the  shock  from  Its  value  ahead 
of  the  shock  (13): 

!♦<*,.  t  .y)2  ■ 

-  I>ox(*;.y)l  /“♦<*;.« .e)de.  (8) 


This  must  be  Integrated  in  conjunction  with  (6a). 
The  ADI  procedure  is  again  adopted  as  outlined  In 
Ref.  (13),  to  effect  a  solution  of  Eq.  (6)  in 
conjunction  with  Eq.  (8),  subject  to  the  time- 
linearized  boundary  conditions  at  Eqs.  (6a)  and  (6b). 

INDICIA!  RESPONSE 


As  a  consequence,  any  boundary  condition  imposed  on 
a  |y|  -  constant  boundary  that  is  less  chan  150 
chord  lengths  away  can  contaminate  the  lndiclal 
response  through  a  reflection  from  a  boundary.  Our 
experience  has  been  chat  an  erroneous  boundary 
condition  such  as  $  *  0  has  to  be  imposed  at  |y| 
greater  than  80  chord  lengths  in  order  to  avoid 
errors  in  the  phase  lag  determined  from  an  lndiclal 
response.  The  same  must  be  true  for  the  computation 
of  a  harmonic  motion,  although  It  would  be  more 
difficult  to  determine  chat  the  phase  lag  was  in 
error  in  such  a  computation.  This  same  observation 
should  also  be  applied  to  unsteady  wind  tunnel 
tests.  If  there  are  significant  acoustic  reflect¬ 
ions  from  the  wind  tunnel  walls,  the  observed  phase 
lags  may  be  in  error.  This  experimental  difficulty 
warrants  further  investigation,  especially  in  two- 
dimensional  studies. 

On  the  other  hand,  we  know  chat  with  the 
appropriate  steady  state  value  of  the  poeenclal 
applied  at  about  twenty  chord  lengths,  the  steady 
state  solution  is  perfectly  adequate.  With  the 
imposition  of  the  unsteady  boundary  condition  (4), 
or  its  time-linearized  analog,  we  find  chat  once 
again  20  chord  lengths  will  suffice.  For  low  to 
moderate  reduced  frequencies,  viz.,  k  ■  0.1  to  1.0, 
the  acoustic  wavelengths  associated  with  the  motion 
are  about  1.0  to  10  chords.  The  grid  spacing 
employed  mey  be  stretched,  but  grid  spacing  compar¬ 
able  to  or  larger  than  the  acoustic  wavelength  will 
result  in  acoustic  reflections  from  the  grid  itself. 
Thus,  while  a  grid  stretching  is  employed  in  the 
calculations,  the  largest  grid  spacing  used  remeins 
a  fraction  of  a  chord  length. 


One  of  the  major  advantages  of  time- lineariz¬ 
ation  is  chat,  for  a  given  mode  of  motion,  the 
amplitude  and  phase  lag  of  the  lift  or  moment 
coefficient  for  a  given  reduced  frequency  may  be 
computed  by  a  linear  superposition  of  the  results 
obtained  for  a  step  change.  For  example,  if  Che 
change  in  lift  coefficient  as  a  function  of  time 
for  a  step  change  In  angle  of  attack,  C ,  ,  is  chat 

3 

sketched  in  Fig.  2,  Chen  the  lift  coefficient  for 
an  angle  of  attack  variation,  a(c),  is 

C,  (C)  -C.  (t)a(o)  +  /*C,  (T)~^  ~  dr.  (9.) 
X  3  a 

Thus,  cor  a  periodic  motion  o(t)  *  a  + 

o 

CtM  -  («)-Ci  (r)2._l“T<i£[  (9t» 

a  3  o 

+  ci  «->v 

3 

The  law  frequency  approximation  made  In  Eq. 

(1),  viz.,  that  k**  was  negligible,  is,  of  course, 

not  valid  for  the  high-frequancy  components  of  tha 
lndiclal  response  calculation.  It  is,  however, 
perfectly  satisfactory  for  the  computation  of  the 
lndiclal  response,  provided  this  response  is  only 
used  to  compute  motions  for  which  k  *  o(l).  In 
order  to  calculate  the  response  for  low  reduced 
frequencies,  however,  we  must  accurately  resolve 
the  lndiclal  response  as  the  motion  approaches  its 
asymptotic  state.  Typically,  this  raqulres  the 
computation  of  the  lndiclal  responsa  for  300  chord 
langths  of  airfoil  motion.  In  this  time  the  un¬ 
steady  perturbations  have  travelled  a  little  more 
than  300  chord  langths  normal  to  the  frecstream. 


Harmonic  Oscillator  Modeling 

Typically,  the  lndiclal  response  of  the  lift 
coefficient  to  a  step  change  In  angle  of  attack, 
flap  angle,  or  imposition  of  plunging  velocity,  is 
like  that  as  shown  in  Fig.  2.  The  same  is  also 
approximately  true  for  the  moment  coefficient  taken 
about  the  airfoil’s  leading  edge  (Fig. 3)  .This  suggests 
that,  to  a  first  approximation,  the  response  is 
nearly  exponential  and  governed  by  a  simple  first- 
order  differencial  equation.  And,  further,  that  to 
a  second  approximation,  the  difference  between  this 
response  and  an  exponential  function  can  be  modeled 
by  a  damped  harmonic  oscillator.  Thus,  if  we  let 
u(t)  be  a  normalized  lndiclal  response  such  that 
u(0)  •  -1,  u(-)  •  0,  we  should  write 

u(t)  -  1  +  u  (t)  ♦  u. ( t )  +  (10) 

0  1 

where 

L  «  =  u  +  Xu  ■  0 

3  0  0  3 

and,  in  general, 

L1«1  5  +  2ptu  +  q^  »  0. 

The  constants  X ar*  determined  to  best 

model  the  lndiclal  response.  Thee  is,  the  solutions 
to  these  equations,  viz., 

u3(t)-.^t  (12a) 

and 

ut(t)  -  ut(o)*’Pitsii»<fle)/a,  (12b) 


3 


where 


a  -/(q1  -  PJ2). 

are  comb in* d  eo  b«sc  approximate  the  indiclal 
response.  To  be  specific,  if  we  let  eQ  *  u(t)  - 
uQ(t),  then  we  choose  X  such  that  J 

y*>  ■  <13) 

is  minimum.  Setting  3 1^/3 X  »  0,  we  find 

X-1  -  2/"{l  -  u(t)}2dt.  (14) 

0 

In  an  analogous  manner  we  let  -  u(t)  -  1  -  u^(t) 
-u^(t)  and  choose  and  q^  so  that 

is  minimized.  This  gives 

p,  -  gw--*?2  and ,  .  - j-oi2.i£.  (is) 

1  2/  (u  -  u  Tdt  1  /  (tt  -  u Jdt 

o  o  o  o 

The  extent  to  which  the  simple  first  approxi¬ 
mation  is  justified  for  selected  examples  is  shown 
in  Figs.  (4)  and  (5).  In  many  instances  an  accept¬ 
able  determination  of  the  phase  lag  requires  the 
second  approximation.  This  will  be  discussed  more 
cully  in  (7).  The  constant  X  and  the  IQ(X)  of  the 

first  approximation  can  be  determined  immediately 
from  Eqs.  (13)  and  (14)  as  the  indiclal  response  is 
being  calculated*  If  I  is  not  sufficiently  small, 
then  the  constants  and  q^  can  be  calculated 

from  Eqs.  (16).  If  the  second  approximation  is  not 
judged  sufficiently  accurate  because  1^  is  not 

acceptably  small,  the  modeling  is  abandoned  as  the 
computational  expense  of  computing  p^  and  q^  Is 

comparable  to  chat  required  for  three  reduced 
frequencies. 

Some  time  ago  it  was  noted  by  Tij deman  (15)  and 
the  authors  .(13)  that  the  amplitude  of  a  harmonic 
response  decays  like  with  increasing  k.  A 

somewhat  more  general  result  is  implied  by  Sqs.  (12). 
If  F(c)  is  a  harmonic  response,  e.g.,  (t),  then 

a 

with 

F(t)  *  F  (t)  +  F  ( t)  +  •••,  (17) 

O  X 

we  find  that  the  first  approximation  gives 

F  /A  •  (1  *  V2)1*  sin(kt  -  9  )  (18a) 

0  0 

where 

slnd  -  k’/(l  +  k,2)1/2.  (18b) 

o 

Here  A  is  the  amplitude  of  the  indiclal  response  to 
a  unit  change  and  k'  -  k/X .  The  second  approxima¬ 
tion  gives 

Fj/A  -  u1(o)k'alsin(kt  -  * 

Ql  +  (k'  +  !5,)2>U  +  (k'  -  3')20l/2  (194) 

where 


slnflj^  -  (k'2  -  a'2  -  1)  i 
(Q.  +  (k’  +  +  (k'  -  a' f~}}112  (19b) 

and  12*  *  ftA  . 

We  see  immediately  from  Eq.  (18)  chat  in  the 
first  approximation  the  amplitude  of  the  harmonic 
response  is 

S  +  (kA^3"1/2 

as  indicated  in  Fig.  4,  which  behaves  like  k  1  for 
large  k,  and  that  the  phase  lag  Is 

sin'1((k/X)/(l  +  (kA)V/2 

which  grows  linearly  with  kA  for  small  k/X ,  and 
thereafter  is  nearly  Independent  of  kA ,  as  can  be 
seen  from  Fig.  (5). 

We  limit  our  discussion  to  the  simple  variation 
of  the  augxlitude  and  phase  lag  of  the  lift  and  moment 
coefficients  for  an  NACA  64A006  in  pitch  with 
reduced  frequency.  Earlier,  more  detailed  results 
depicting  the  3hock  motion,  etc.,  are  to  be  found 
in  Refs.  (13)  and  (16).  Because  we  have  linearized 
about  a  steady,  small  perturbation  solution,  we 
draw  no  practical  conclusions  from  our  study. 

Tij deman  (personal  communication)  reports  that  the 
application  of  LTRAN2  to  determine  the  response 
about  an  experimental  steady  state  for  the  F29  air¬ 
foil  at  varying  incidence,  in  conjunction  with  strip 
theory  and  the  results  of  panel  methods  for  suberi- 
tical  flow  to  account  for  three-dimensional  effects, 
was  successful  in  predicting  the  flutter  boundary 
of  the  F29  wing. 

Figures  4  and  5  give  the  amplitudes  and  phase 
lags  of  the  lift  and  moment  coefficients  for  an 
HACA  64A006  airfoil  oscillating  in  pitch  at  selected 
reduced  frequencies  with  *  0.86  and  0.88.  Indi¬ 
vidual  results  are  shown  by  symbols  with  the  reduced 
frequency  noted  below  them.  Generally,  they  are 
well  described  by  the  first  approximation  of  the 
harmonic  oscillater  model.  For  -  0.36  the  lift 
and  mid-chord  moment  results  are  indlstlnqulshable, 
but  their  phase  lags  are  not  correctly  captured  by 
the  modeling  of  the  first  approximation.  In  this 
case,  the  Initial  part  of  the  indiclal  response  Is 
not  correctly  captured.  This  is  easily  dealt  with 
in  the  model  without  going  to  the  second  approxi¬ 
mation  (7). 

COMPUTATIOKAL  EFFORT 

As  noted  earlier,  a  time- linearized  calculation 
requires  substantially  less  computational  effort 
chan  a  nonlinear  one.  We  delineate  those  differ¬ 
ences  here,  noting  that  the  larger  grid  spacing  that 
may  be  used  with  our  shock-fitting  procedure 
Implies  a  computational  saving  in  addition  to  the 
considerations  discussed  here. 

If  we  let  T  represent  the  number  of  time  steps 
required  to  calculate  a  time  periodic  solution  In 
an  l*M*N  special  domain,  then  the  total  computa¬ 
tional  effort  for  a  flutter  study  using  either  a 
nonlinear  or  time- linearized  algorithm  Is  propor¬ 
tional  to  the  product  TLMN.  If  we  do  not  time- 
linear  ize,  this  must  be  done  at  K  reduced  frequen¬ 
cies  to  give  a  total  computational  effort  that  is 
proportional  to  KTLMN.  If  a  time- linearized 
procedure  Is  used  to  compute  a  single  indiclal 


response,  Che  effort  Is  TLMN.  This  can  be  used  to 
generate  the  mode  response  for  any  reduced  frequen¬ 
cy  in  about  T2  additional  steps,  giving  a  total 
computational  effort  proportional  to  TLMN  +> 

(const) *T2K.  Typically,  T  is  500  and  L,M,N  are 
about  50,  25,  25,  respectively.  Thus,  a  nonlinear 
analysis  at  ten  reduced  frequencies  has  a  computa¬ 
tional  effort  of  about 

10- 500- 50- 23- 25  »  108. 

On  the  ocher  hand,  a  time- linearized  computation 
requires  a  computational  effort  of 

500*50.25*25  + 

(const)  (500)210  *  I07  +  (conse)'l 06 

providing  a  factor  of  K  reduction  over  the  nonlinear 
analysis.  The  constant  of  proportionality,  and  the 
less  refined  L  grid  spacing  required,  also  favor 
the  time- linearized  algorithm,  contributing  roughly 
another  factor  of  ten  reduction  over  the  nonlinear 
procedure. 

In  two  dimensions,  the  two  computational  compo¬ 
nents  of  the  time- linearized  calculation,  viz.,  the 
indlcial  response  and  Its  linear  superposition  for 
each  reduced  frequency,  are  comparable  and  an 
efficient  integration  algorithm  must  be  used  for 
Eq.  (9b). 

CONCLUSION 

The  amplitudes  and  phase  lags  of  the  lift  and 
moment  coefficients,  at  selected  reduced  frequencies, 
can  be  computed  accurately  and  efficiently  by  cime- 
llnearlzaeion  about  a  measured  or  computed  pressure 
field.  For  three-dimensional  studies  the  linear 
superposition  of  the  results  of  a  single  indlcial 
response  substantially  reduces  the  computational 
effort.  Further  reductions  are  achieved  by  using 
the  linearized  far-fleld  for  the  unsteady  flow,  and 
by  creating  the  shock  wave  as  a  discontinuity  in 
Che  computations.  In  many  cases,  the  indlcial 
response  can  be  modeled  by  a  simple  harmonic  oscil¬ 
lator,  and  this  provides  an  analytical  result  for 
the  dependence  of  the  lift  and  moment  coefficients 
on  reduced  frequency,  further  reducing  the  compu¬ 
tational  effort. 
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TIME  in  chords  travelled 


Figure  3.  Moment  coefficient  as  a  function  of  time  for  a 

step  change  in  angle  of  attack;  NACA  64AQ06,  *  0.88. 
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Figure  4.  Lift  and  moment  amplitude*  normalized  by  their  quaai-steedy  value  as 
function  reduced  frequency  for  an  NACA  64A0Q6  airfoil. 


